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Abstract 

A new approach to finite basis sets for the Dirac equation is developed. It 
solves the problem of spurious states and, as a result, improves the conver- 
gence properties of basis set calculations. The efficiency of the method is 
demonstrated for finite basis sets constructed from B splines by calculating 
the one-loop self-energy correction for a hydrogenlike ion. 
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At present, a great variety of calculations in atomic physics and quantum chemistry 
are based on finite basis sets. First successful attempts to utilize finite basis sets in rel- 
ativistic quantum mechanics were made many years ago [1-6]. Application of B splines 
for constructing basis sets [7-9] provided new impact to this field. Nowdays, B sphnes are 
widely employed in computational atomic and molecular physics (see [10,11] and references 
therein) . 

In contrast to the nonrelativistic case, the use of finite basis sets in the relativistic 
theory is generally accompanied by the occurrence of spurious states (see, e.g., [4,9]). For 
the attractive Coulomb potential, spurious states appear for > as the lowest bound states 
with non-physical energies {k — {—iy~^^~^^^^{j + 1/2) is the quantum number determined by 
the angular momentum and the parity of the state). For the point nucleus, their energies 
coincide with the physical energies of the lowest bound states with k. < 0. The wave functions 
of these states oscillate rapidly and, therefore, in many cases they may be disregarded in 
practical atomic calculations. However, since the presence of the spurious states disturbs the 
spectrum, it worsens the convergence properties of the basis set calculations in some cases. 
For this reason, despite of the early promising results [13,14], the finite basis set method 
has not been extensively employed in calculations of radiative corrections. To date, most of 
these calculations are performed by means of analytical or numerical representations for the 
Coulomb-Green function (see, e.g., [15-18]) or by the space discretization method [19-21], 
in which the spurious states are eliminated from the very beginning (see [19] for details). 
The problem of spurious states is especially demanding in calculations of atoms in strong 
external fields and for molecules, where, generally speaking, they can not be selected as the 
lowest bound-state energies and, therefore, can not be eliminated without employing special 
methods. Furthermore, even in those cases where spurious states do not play any significant 
role, their presence generates difficulties for providing an adequate estimate of the accuracy 
of the calculations. 

A number of proposals for solving the problem of spurious states were presented previ- 
ously [22,9,23-25]. Our analysis of these methods indicates, however, that none of them can 
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be considered as completely satisfactory. We find that treating the problem by a particular 
choice of the boundary conditions in the B-spline method [9] does not provide any improve- 
ment in calculations of radiative corrections. The other methods either require considerable 
modifications of the standard numerical procedure or they are strongly limited to a specific 
choice of finite basis sets. Their apphcability to calculations of radiative corrections thus 
remains questionable. For this reason, we suggest a new method for solving the problem of 
spurious states. The efficiency of the method is demonstrated for finite basis sets constructed 
from B splines by calculating the one-loop self-energy correction for a hydrogenlike ion. 
For the case of a central field V{r), the Dirac wave function is conveniently represented 

by 



y iF{r) n-^m{T^) 

where n = r/r. With this representation, the radial Dirac equations can be written as 



(1) 



(2) 



where (in units: ^ = 1) 



mc +V c 
d ^ K 
dr r 



and 



d K-i \ 

dr r 



—mc^ + V 



(3) 



0(r) 



^ G{r) ^ 



(4) 



is the two-component radial wave function. The scalar product of the two-component func- 
tions is defined by 



(0a|0fe) = / dr [Ga{r)G,{r) + F,(r)F,(r)] . 
Jo 



(5) 



The radial Dirac equations can be derived from an action principle 5S — with 
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S ^ {<t>\H,\<t>) - E{4>\<t>) , 



(6) 



if proper boundary conditions for G{r) and F{r) are implemented. The functions 0(r) can 
be approximated by 



2n 



(7) 



i=l 



where the two-component functions rtj(r) are assumed to be square integrable, hnearly in- 
dependent, and satisfying proper boundary condition at r = 0. The variational principle 
reduces to the following algebraic equations 



dS/dci^O, i = 1,2, ...,2n. 



(8) 



This leads to a generahzed eigenvalue problem 



KikCk — EBikCk , 



(9) 



where — {{ui\Hi^\uk) + {uk\Hi^\ui))/2, — {ui\uk), and the summation over repeated 
indices is implied. 

Let us first demonstrate that the widely applied choice 



Ui{r) = 



Ui{r) 



, i = 1, ...,n, 



V » / 

^ ^ 



, i — n + l, 2n , 



(10) 



(11) 



where {7ri(r)}"^i are square integrable functions satisfying the boundary condition 7ri(0) = 0, 
results in the occurrence of spurious states. In this case, equation (9) takes the form 



(mc^ + V - E)ikPk + cDikQk = , 
c{D%pk + {-mc^ + V- E)ikqk = , 



(12) 
(13) 



where (±mc^ + V — E), D, and D'^ are nxn matrices with elements 
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{±mc'' + V-E)ik^ / dr ni{r){±mc^ + V - E)7rk{r) , (14) 

Jo 

poo / d K\ 

Dik^ dr7ri{r)[-— + -)7rk{r), (15) 
{D% = dr7r,(r)(^ + -)7rfe(r), (16) 

and Pi = Ci, qi = CiJ^n for i = 1,2, ...,n. Let us consider the nonrelativistic limit (c — > oo) 
and introduce vectors P and Q with components {pi}"^i and {?i}"=i, respectively. Then Eq. 
(13) yields Q — {l/2mc)D'^P. Substituting this expression into Eq. (12), we obtain 

DD^P + 2m(mc^ + V - E)P = . (17) 

For the pure Coulomb field, V{r) — —X/r (A > 0), introducing the matrix Cik — Dik — 

{m\/ K)5ik, Eq. (17) reduces to 

C,ClP = eP, (18) 

where e = \2m{E — mc^) + rn^X^ / k^] and the dependence of the C matrix on k is exphcitly 
indicated. Taking into account that C-^ = we find that the corresponding equation 

for k' — —K can be written as 

ClCP'^e'P'. (19) 

Multiplying Eq. (18) with C^, we obtain that each nonzero eigenvalue of C^C^ is an eigen- 
value of C^C^. Evidently, the inverse statement can be proven in a similar manner and the 
dimension of a nonzero eigenvalue subspace is the same for C^C}. and C}.C^. Accordingly, 
the spectra of C^C^ ^I^k may differ only by the dimension of the zero eigenvalue 
subspace. For finite matrices, the dimension of the subspace with e = is the same for 
CkC\ and ClC^, since the total number of eigenvectors as well as the dimension of the 
nonzero eigenvalues subspace is the same for C^C^ and C^C^. Therefore, the finite matrices 
CkC}. and C^C^ = C-^C^-k have an identical spectrum. Conversely, we know that the exact 
analytical solution of the Dirac equation for the Coulomb potential yields different lowest 
bound-state energies for k < and ac > 0. This is due to the fact that within the exact 
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(infinite dimension) treatment the subspace with e = may have different dimensions for 
K, <0 and K, > cases. This can easily be checked by solving the equation 



which in case of finite dimensions is equivalent to the equation C^P — 0. Solving equation 
(20) yields G{r) = Aqt''^ exp {mX/ nr). For «; < 0, this solution has the proper behaviour 
at r — > and at r — > oo. However, this does not hold for k > 0. Concluding, in the 
approximation of finite dimensions, our proof clearly indicates the presence of spurious 
states with energy E — mc^ = — mA^/2«;^ corresponding to e = 0. It can be shown that this 
result remains valid for the full relativistic theory as well. It is obvious that spurious states 
must occur for any other potential one is dealing with in atomic calculations. 

To eliminate the spurious states, in Ref. [23] "kinetically balanced" Slater type functions 
were employed. Within this method, for k > the lower components in equation (11) 
are replaced by functions Pi{r) which, in the nonrelativistic limit, are related to the upper 
components 7rj(r) in equation (10) via 



This method provides a high accuracy in calculations of bound-state energies in atoms for 
both sign of k, if an extended nuclear charge distribution is introduced and proper bound- 
ary conditions are implemented [23]. However, since the basis functions are "kinetically 
balanced" for positive energy states only, the application of this method to calculations of 
the QED corrections might be problematic [26]. The equivalent treatment of positive and 
negative energy states would provide, in particular, the well-known symmetry properties of 
the Dirac spectrum under the transformations V — >■ —V, k —k,, and G ^ F. It is evident 
that the "kinetically balanced" functions do not meet with this requirement. 

In the original version of the B-sphne method [8-10], to achieve that the first positive- 
energy states K, > correspond to physical bound states, an additional term had to be 
introduced in the Hamiltonian, which formally implements the so-called MIT boundary 



{d/dr + k/t - m\/K)G{r) = , 



(20) 



Pi{r) ~ {l/2mc){d/dr + K,/r)7Ti{r) . 



(21) 
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condition [12]: G{R) — F{R), where R is the cavity radius, together with the condition 
G{0) = 0. However, since the presence of the additional term does not imply any practical 
advantages, it is usually omitted in calculations. Instead, the boundary conditions are 
generally implemented by eliminating the first and the last basis function which are the only 
ones that do not vanish at r = and r — R, respectively. This method was successfully 
employed for calculations of the two-photon exchange diagrams within the rigorous QED 
approach [27-31] and for relativistic calculations of the recoil effect [32-34]. However, its 
application to calculations of pure radiative corrections [13,14] was less successful, compared 
to the other methods [15-18,20,21]. We conjecture that this would not be the case if the 
spurious states were eliminated in a more natural manner than it was done in [8,9]. 

It is known (see, e.g., [14,10,23]) that the case of the pure Coulomb potential requires 
generally a special care in implementing finite basis set methods. This is due to the singu- 
larity of the Coulomb potential at r — > 0. However, in practical calculations it is standard 
to modify the potential to account for the finite nuclear size, which eliminates this problem. 
For this reason and for simplicity, we restrict our consideration to the case of a finite nu- 
clear charge distribution, bearing in mind that the limit of a point nucleus can be treated 
by extrapolating a series of calculations for extended nuclei to vanishing nuclear size. For 
extended nuclei we propose to employ the following basis set 



where the linearly independent functions {7rj(r)}"^i are assumed to be square integrable and 
satisfying the proper boundary condition at r = 0. We state that this basis set satisfies the 
following requirements: 

1. It is symmetric with respect to the replacement k ^ —k and the interchange of the 
upper and lower components. 




(22) 



(23) 



V 



7 



2. The functions ui, Un provide the correct relation between upper and lower compo- 
nents for \E — mc^l, |V^(r)| <^ 2mc^, while the functions Un+i, U2n do the same for 
\E + 7710^1, \V{r)\ < 2mc^. 

3. Calculations utilizing the standard finite basis set determined by Eqs. (10) and (11) 
can be easily adopted when employing the basis (22)- (23). 

4. No spurious states occur for attractive as well as for repulsive potentials. 

The properties 1-3 follow immediately from definitions (22) and (23). The absence of 
spurious states can be explained as follows. Performing similar steps as for the derivation 
of Eq. (17), for \E — mc^\ ^ 2mc^ we obtain 

— LP + (y + mc^ - E)P = , (24) 
2m 

where 

= fdr <r){-^^ + + ^H') ■ (25) 

Eq. (24) takes the form of the ordinary Schrodinger equation with / = |«;| -|- {k,/\k,\ — l)/2 
in the finite basis representation. As is known, it generates no spurious states. The region 
\E + mc^\ <^ 2mc^, where spurious states may exist for repulsive potentials and for k < 0, 
can be considered similarly. In this case, we obtain the equation 

MQ + {-V + mc' + £;)g = , (26) 

2m 

where 

M. = f*..M(A + l)(-^ + ^).,W. (27) 

Eq. (26) has also the form of the ordinary Schrodinger equation but with V — \k\ — (k/|k| -I- 
l)/2. It transforms into equation (24) under the replacements n — >• — V —V, E — —E, 
Q ^ P and does not generate any spurious states. This is a consequence of the equivalent 
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treatment of the positive and negative energy states. For this reason, the new basis may be 
termed conventionally as dual kinetic-balance (DKB) basis. 

The validity of statement 4 has also been proven by numerical calculations with 7ri(r) — 
Bi{r), where -Bj(r) are the B splines defined on the interval (0, i?) as in Ref. [9]. The first 
and the last spline function have been omitted. Standard test calculations (see, e.g., [9]) 
show that this basis satisfies suitable completeness criteria as n — > oo. 

Finally, let us consider the calculation of the one-loop self-energy (SE) correction to 
the ground-state energy of a hydrogenlike ion employing the new basis set. Generally, the 
SE correction is expanded into the zero-, one-, and many-potential terms. The ultraviolet 
divergences in the zero- and one-potential terms and in the counterterm cancel each other 
and their evaluation can be performed according to the formulas presented in Refs. [35,36]. 
As to the many-potential term, although it does not contain any ultraviolet divergences, its 
calculation is most difficult since it involves the summation over the whole Dirac-Coulomb 
spectrum. In Table 1, we compare our results obtained for the many-potential term for 
Z = 20 employing the DKB basis set (22), (23) with 7ri(r) = Bi{r), the old basis (10), (11) 
with the same 7ri(r), and the results of a calculation using the analytical representation for 
the Coulomb-Green function. The shell model for the nuclear charge distribution has been 
used with the radius R = 3.478 fm. In the basis set calculations, the contributions with 
\k\ > 10 were obtained by an extrapolation. Adding the zero- and one-potential terms to the 
many-potential term yields 0.06409 a.u. for the old basis, 0.06426 a.u. for the DKB basis, 
and 0.06425(1) a.u. for the Coulomb-Green function calculation. This comparison clearly 
demonstrates a significant improvement in accuracy, if the DKB basis is employed instead 
of the old one. More extensive calculations employing the DKB basis will be presented in 
forthcoming papers. 
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TABLES 

TABLE I. The partial-wave contributions to the many-potential term (in a.u.) for the Is state 
at Z = 20, obtained by the basis set methods and by the Coulomb-Green function (CGF) method. 
The number of the basis functions: n = 60. The shell model for the nuclear charge distribution 
has been used with R = 3.478 fm. 



Term 


Old basis 


DKB basis 


CGF method 


\k\ = 1 


0.848691 


0.848750 


0.848741 


\k\ = 2 


0.020618 


0.020662 


0.020653 


K = 3 


0.005302 


0.005331 


0.005326 


\k\ = 4 


0.002121 


0.002139 


0.002137 


|«;| = 5 


0.001050 


0.001062 


0.001062 


\k\ = 6 


0.000590 


0.000597 


0.000598 


|«;| < 9 


0.879127 


0.879303 


0.879288 


\k\ > 10 


0.000587 


0.000585 


0.000583 


Tot. many-pot. 


0.87971 


0.87989 


0.87987(1) 
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